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Abstract 

Solvable via Bethe Ansatz (BA) anisotropic statistical model on cubic lattice consisting of locally inter- 
acting 6-vertex planes, is studied. Symmetries of BA lead to infinite hierarchy of possible phases, which is 
further restricted by numerical simulations. The model is solved for arbitrary value of the interlayer coupling 
constant. Resulting is the phase diagram in general 3-parameter space. Two new phases of chiral (spiral) 
character and new first order phase transition appear due to the interplane interaction. Exact mapping onto 
the models with some inhomogenious sets of interlayer coupling constants is established. 



1 Introduction 

Exactly solvable models e.g. models for which set of physical quantities such as the bulk free 
energy, the interfacial tension, some critical exponents etc. can be calculated analytically, 
play important role in statistical mechanics and the phase transitions theory. Most solvable 
statistical models are two-dimensional (see for review). YBE, or star-triangle equation, 
serves as integrability condition, or transfer-matrices commutativity condition for 2D solvable 
models. Unlike 2D models, solvable models in are very rare examples. Usually solvable 3Z) 
models — see e.g. Zamolodchikov model and its generalizations 0] — are based on solution 
of tetrahedron equation. The latter is natural generalization of Yang-Baxter equation YBE and 
serves as transfer-matrices in 3D commutativity condition. 

Unfortunately all 3D solvable models known so far possess common rather unsatisfactory 
feature, from the viewpoint of applications to statistical mechanics — they incurably have 
negative or even complex Boltzmann weights. Since these weights are probabilities (up 
to a normalization), this property makes statistical mechanical interpretation of the models 
highly problematical (although the associated quantum problem may still be sensible). 

Recently the method was proposed in |p, ^ which allows to construct solvable statistical 
models (however anisotropic ) with positive Boltzmann weights and local interactions in 3D 
starting from solvable 2D models. The idea is the following: the whole row of infinite number of 
separate sites (vertexes) is considered as a simplest object. The states probabilities (Boltzmann 
weights) are defined as product of those for each site separately, multiphed by interaction factors 
depending upon the local configurations of each two neighbouring sites. Such construction leads 
to multilayered models on 3D-lattice, consisting of 2D integrable planes-layers, with specific 
interaction between them. What is important is that transfer-matrix commutativity condition 
turns out to be not tetrahedron equation but an infinite set of usual Yang-Baxter equations. 
These equations can be satisfied simultaneously, producing the variety of 3D-extended mod- 
els satisfying the necessary physical requirements — positivity of Boltzmann weights and 
locality of interactions. 

Here we investigate the simplest example of such a model — 3D-extended 6-vertex model, 
with the homogenious set of interplane coupling constants. 6-vertex model is rather popular 
object to study, because of being limiting case for many solvable models in critical region (see 
e.g. [jl|). 3D solvable extension for 6-vertex model was obtained in and in general- 
ized onto other solvable vertex models with charge conservation and inhomogenious sequences 
of interplane interaction constants. In the phase diagram for the free fermions case was 
conjectured. Here we go beyond this restriction and obtain the phase diagram in the general 
3-parameter space, for the extended 6-vertex model. 

The paper is organized as follows: first we remind the definition of the anisotropic 3D- 
extended 6-vertex model and consider the strong interplane interaction limit. Then we give the 
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Figure 1: Multivertex consisting of K simple vertexes 



alternative derivation of Bethe Ansatz (BA) equations, using the established gauge equivalence 
of our multilayer model to the set of 2D 6-vertex planes, each one in a field defined by the 
polarization of planes - neighbours. Using the symmetries of BA, we prove the equivalence 
between the model with homogenious set of interaction constants and some models with inho- 
mogenious sets. We make use of these results, together with hypothesis of non-degeneracy of 
maximal eigenvalue, to eliminate the problem. The resulting phase diagram is obtained in the 
next section. Conclusion and discussing the possible generalizations closes the paper. 



2 Definition of the model 

The model we consider is system of K planes. Each of these planes is the symmetric 6-vertex 
solvable model on square N ^ M lattice ^j. We can "paste" together ij-sites of all planes and 
formally get 2D system with complex site consisting of K simple vertexes (see fig. 1). 

Boltzmann weight of such multivertex fragment has form 

L:;:^;;"J§g Jf = n L.^'^jfie^vi-huaA^, + h,a',^M (1) 

k=l 

where hk are arbitrary constants, defining interaction between nearest neighbours in k- 
th and {k + l)-th plane, Lg^^^'f? — Boltzmann weights of "source" 6-vertex solvable model, 

state variables ak,f3k, ■ ■ ■ take values ±1. For the planes (layers) we impose periodic boundary 
conditions K + 1 = 1. The permissible configurations of 6-vertex model are drawn on fig.2. The 
variables sitting on edges are represented by arrows on fig.2. "-I-" ("-") correspond to arrows 
pointing up or to the right (down or to the left). The local Boltzmann weights are invariant 
under inversion of all arrows. 



Strong interplane interaction limit h ^ 00. 

To understand the nature of interaction for the homogenious model (all interplane interaction 
constants hk equal: hk = h, any k), let us consider the strong interplane interaction limit 
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Figure 2: Fig. 2. Permissible vertex configurations for 6-vertex model. Boltzmann weights of configurations: 
wi = W2 = a; W3 = W4 = b; = wq = c 

h oo. The configuration with maximal Bolzmann weight will give the main contribution to 
the partition sum. To obtain this configuration, we should maximize exponential factor in (|l|) 
for all sites in all planes. 

Without loss of generality, we can set = = +1 for A;-th plane. Then, to maximize 
the exponential factor in (|l]), we should choose Pk+i = —1; (^k+i = +1 in + 1-th plane, 
i.e. vertex of type 4 (see fig. 2), and a^+i = +1; = —1. Performing the next step, we 

get Pk+2 = ftfc+i = —1- This fits vertexes of type 5 and 2; however vertex 5 is unsuitable^ 
Proceeding analogously further, we get: planes k — {k + 4) are formed by the vertexes of type 
1, 4, 2, 3, 1, respectively. It is easy to see that each next vertex (next plane) can be obtained from 
the previous one by the clockwise rotation by 7r/2. After full rotation over 27r, the configuration 
repeats. 

Thus, in the strong interplane interaction limit, the model has homogenious structure within 
each horizontal plane, and spiral structure with period 4 in vertical direction: each next plane 
configuration is obtained from the previous one by clockwise 7t/2 rotation: 

. . . / \ / \ . . . (2) 

The arrows in formula above are the resulting from vector summation of all arrows in each 
plane: for example, the arrow y corresponds to homogenious configuration of vertexes of type 
1 in Fig.2. 



3 The matrix formulation 



It is convenient to rewrite (|l|) in matrix form. Let us introduce 2^^ x 2^^ matrix £^ acting in 
the tensor product Oi^ '^fl'o Hi ®9n with elements 

/ 0.-10.2 ■■■Cii.K 0102 •••0K rrrr Si £ { 'L iQ \ 'L f Qf \ 

\^ Ja;a^...«'^/3;/3^.../3^ = ^ a^a\0 p^fi'^d a^a'^Q fi^fi'^ ■ ■ ■ Oak_^a'^^__^0 f3^_^f3'^__^ CXp [-llkakPk+1 + 

Lq^ is the 4x4 local L-matrix of six-vertex model with the following nonzero elements: 

-^11 ~ -^22 ~ '^j -^12 ~ -^21 ~ -^21 ~ -^12 ~ ^ (4) 

in which a, b, c are Boltzmann weights of permissible configurations (see fig.2). Here and below 
in this section we borrow notations from paper of Faddeev and Takhtajan ||^. Note that 

^ In order to get maximal Boltzmann weight, we have to construct configuration in each plane with the same type of vertexes. 
This can be done with type-1 — 4 vertexes but not with type-5 vertexes -they always go in pairs with type 6 vertexes. 
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matrices ^ corresponding to neighbouring /c-s, do not commute. The multivertex Boltzmann 
weights (|l]) are matrix elements of the ordered matrix product of over all planes: 

K 

fc=l 

K 

= H -^6t;a^;^^exp{-/ifcafc/5fc+i + } (5) 

k=\ 

(expressions (|I]) and (|^) are the same, due to periodicity). 

Using introduced notations we can write partition function of the model in usual (see e.g. 
0) form: 

Z = Tr(T)^ (6) 

with transfer-matrix T being the trace of the ordered product of local matrices L along the 
row (we'll denote matrix L acting in the site n as L„): 

N 

T = rr(nL„) (7) 

n=l 

Here Tr operation and matrix product goes over a-indexes (we'll call it: 'auxiliary' space) and 
Tr operation in the previous formula (6) goes over /5-indexes ('quantum' space). 

The free energy per site in the thermodynamic limit iV, M, K ^ oo'vs, defined by the maximal 
transfer-matrix eigenvalue 

T^ = M>- f = -kBT \nn -^\A^ma.\. (8) 

Important feature of our i^'-plane model is that its monodromy matrix T can be written as 
an ordered product of more simple ones (we'll denote matrix from (^ acting in the site n 
as £^): 

N N K K N K 

^ = n Ln = 11(11 ^n) = n(n ^n) = n (9) 

n=l n=l k=l k=l n=l k=l 

N 

Tk=l[^n (10) 

n=l 

Matrix C'^ from (^) can be written in compact form as can be easily verified: 

C'' = exp (-/ifca('^V('=+i))L6/ exp (/i^.^r^'^V^'^-i)) (11) 

where we denote by the o"'-^^ and r'^'^^ the diagonal matrix = diag(l, —1) acting nontrivially 
in k-th 'auxiliary' and k-th 'quantum' space, respectively: 

/ -.(k)\oiia2...aKl3il32---I^K _ £ £ ( f f y^K r 

Ja',a^...a^/3i/3^.../3^ - °o,^a'^ ■ ■ ■ <^ak-i»',^_S^ ) a^a'^f^ a^+ia'^^^ • • • "^ok^ lli=l 

(12) 

( ^(k)\oi\Oi2...aKPiP2...PK _ f ^ ( „z\ r r tt-'^ A 



5 



4 Bethe Ansatz equations 



It can be shown - see Appendix for a proof - that the 3D model under consideration (|l|) is 
gauge equivalent to the set of 2D 6-vertex planes, each in its own horizontal field with the 
strength defined by the vertical polarization in neighbouring planes and interplane interaction 
constants: 

for the fc-th plane the field strength 

Hk = hkVk+i - hk-iVk-i (13) 

where polarization is defined as usual: 

nl-ni 2nl-N 
= ^ = ^ (14) 

n\ (n^) is the number of upward (downward) pointing arrows in horizontal row in plane k. 

Therefore, the Bethe Ansatz of the model is given by the well-known formulas for 6-vertex 
model in an external horizontal field (see e.g. 0); below Uk = 4, yk = ^"^^ ): 

A.(i/.) ^ a-e-"^ fl " - ^' . ^^e-^"^ li " ' "'^^ (15) 

j=i a — bTj j=i —a + 6rj 

(k) 

where r satisfy the set of Bethe Ansatz equations 



nfe Ak) (k) _ c^A ik) , . 

_ r r 2Ar + ^^^^ 

1=1 rj 'tI ' - 2ATj; + 1 



where 



nk = 1,2,... N; k = l,2...K 
+ - 



(17) 



2ab 

a,b,c are Boltzmann weights of symmetrical 6-vertex model configurations (see fig. 2). 
The global transfer-matrix eigen-value is the product of those over all planes: 

= Kl t = KmKm ■ ■ ■ AyAHK); y, = (18) 

The Eqs. ( |r5|JT6|) were obtained directly in using the quantum inverse scattering 
method, and the analytic ansatz method, respectively. 

We shall restrict ourselves to the model with all interplane interaction constants equal: 
hk = h, any k (below we refer to that case as to the homogenious model). However as is 
shown below, our results are valid also for some models with inhomogenious sequences {hk}. 

To obtain the bulk free energy and complete phase diagram, one should find such a sequence 
of Uk, {nk}f^i that corresponds to the maximal transfer-matrix eigen-value Amax, and A^ax 
itself, in complete 3-parameter space a/c, b/c, h. This program can be performed easily in the 
following limiting cases: 

^We thank Yury Stroganov for drawing our attention to this fact 
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1) for the single 6- vertex model (decoupled planes limit h = 0). Then BA solutions t^''^ lie 

on a unit circle rj'^'' = e*'^j , A^'^^ real. The expression for A^ax is found then by the integral 
equation method (see e.g. |1|]). 

2) 'Ising chain' case A = ±oo, i.e. a = or 6 = 0, any h is equivalent to A = ±oo, h =0. 
Thus it is reduced to the previous case. 

3) for the free fermions case A = 0. Then RHS of (|16D is equal to 1, and BA is solved 
trivially. The complete analysis is done in [^]. 

In the general case i.e. nonzero h and A, the structure of BA equations doesn't permit 
simple analysis. The reason is in that case the locus of BA roots is unknown. The problem is 
similar to that arising when one considers usual 6-vertex model in external horisontal field (see 

ii)- 

Connections between models with interplane constants 
{. . . h, h, h, h, ...},{.. . h, -h, h, -h, . . . }, and arbitrary sets of 'h' and '-h'. 

The models listed in the title look different, and different they are, having for instance, 
different 'strong interplane interaction limit' h ^ oo (this can be verified directly as it was 
done for homogenious model at the beginning of the paper). We'll show however that they 
have precisely the same transfer-matrix spectrum, and therefore the same phase diagram, with 
minor redefinition of phases. 

Let us take the state of the fc-th plane with arrows pointing up, in a field Hk, described 
by Eqs. (^,|T^). Reversing of all horizontal and vertical arrows, together with changing the 
sign of Hk, leaves the Boltzmann weights and therefore the eigen- value (|15|) invariant: 

Ay.iHk) = A_y^i-Hk) (19) 

(Bethe Ansatz (p!6D changes accordingly). Using this formula, write the global transfer-matrix 
eigenvalue (|18|) for the set {yi,y2, -ys, -y4,y5,y6, -Vr, -Vs, ■■■ - Vk} for the homogenious 
model: 

A^, ly.^'.y, = AyMy2 + y k)) Ay^'V^ - yi))A-yM-yi " ^2 ) ) A ( /i (^/s + y-i)) = 

Ay,{h{y2 + yK))Ay,{h{-y^ - yi))Ay,{h{yi + y2))AyM-y5 - Vs)) (20) 

Remark. Here and below we assume the number of planes K to be infinitely large and 
divisible by all numbers ii' = 2*3*4*...,i^+l = l, to avoid complications connected with 
the boundary effects. For the sake of simplicity we write down only the significant part of the 
multiplication (|T^), then it continues periodically. 

Let us write down the global eigenvalue for the system with alternating constants: 

K mv'y.'' = KiKy2 + yK))AyM-y^ - yi))KMy4 + y2))AyM-y^ - ys)) (21) 

Comparing with the previous formula we have 

\h h h h \h —h h —h (')')\ 

yi y2 -ys -y4 ~ yi 1/2 ys y4, \ ' 

(periodical continuation is implied - see Remark above). 
Analogously one obtains 

Ah h h h \h h —h —h (')'^\ 

y\ J/2 ys J/4 ~ yi j/2 ys -yi \ ' 
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Actually, one could coin such transformations between the model with homogenious sequences 
and those with arbitrary sequences of 'h' and ' — h': {/ifcjfLj^, hk = ±1, k = 1,2, . . . K. For 
instance, 

\h h —h \h h h h (0A\ 

yi y2 ys ~ yi 1/2 ys -yi -ys -ye ^ ) 

and so on. 



5 Symmetries 

Here we make use of the results of previous paragraph to establish transformations which map 
the eigenvalue to the eigenvalue of the same model. 

First of all, note that simultaneous inversion of polarization in all planes Z 

z{y} = {-y} (25) 

leaves the eigenvalue invariant: 

k k 

This is also the direct consequence of the fact that local Boltzmann weight (|I]) is invariant 
under transformation all ak,Pk —ctk, —Pk- 

Another transformation we obtain, inverting the order of {?/}-set and {—h}-set: 

yi y2— vk yi yx yK-i— y2 

(this is proved analogously) Transformation J 

J{yiy2--- yx} = {yi yx yx-i • • • 2/3 2/2} (26) 

is a symmetry for the model with alternating constants: 

\h —h h —h \h —h h —h 

^\y} - ^^J{y} 

Then, one more independent symmetry exist for that model: 

F{y} = {y2 -y, ... yx -2/1}; ' = A^^-^f ' (27) 

For completeness, we define operator of one-step shifting S: 

S{y} = {2/2 2/3 - • • yK yi} (28) 

Evidently, S*", integer n, is a trivial symmetry for the homogenious model and 5*^" is a trivial 
symmetry for the model with alternating constants. Transformations Z, J, F, S"^ are basic sym- 
metries for the model with alternating constants. Let us obtain the basic symmetries for the 
homogenious case. Z and 5* are the symmetries already. 

Then, define the transformation from (|2^) between the two models: 

<5{2/} = {2/1 2/2 -2/3 -2/4- • •}; Q~^=Q 

Transformations QJQ and QFQ are the symmetries sought for: 

P{y} = QJQ{y} = {2/1 -yK yx-i -yK-2 ... 2/3 -2/2}; = I (29) 
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QFQ = S 



So, for the homogenious model all = h two nontrivial symmetries Z and P exist defined 
by the (^51) , (p9D . Working out the same procedure for the model h, h, —h, —h . . . doesn't give 
additional information. It is noteworthy also that the Z-symmetry is not independent but can 
be expressed in terms of shifting and P-symmetry 

Z = {spy 

Nevertheless it proves convenient to keep it in mind as a separate symmetry. As is shown 
below, Z and P symmetries, having been combined with the non-degeneracy hypothesis, put 
drastic constraints on the polarization set {y} which corresponds to the maximal eigenvalue of 
the global transfer-matrix of the homogenious model. 

6 Finding the maximal eigenvalue set {y} 

We are now in a position to solve the maximal eigenvalue problem, for the homogenious model. 
First we have to point out which set {y} corresponds to it for any value of parameters a/c,b/ c, h. 
We shall parametrize the possible sets by a period T and define the set {y^^^ as that with 
periodically repeating entries: 

{2/^^^} = {yi 2/2 • • • 2/T 2/1 2/2 • • • 2/T • • • Vt}; {y^^^lr+n = {2/^^^}n 

Note that cases T = 1 and T = 2 are trivial because they lead to vanishing of interaction h 
from all expressions (see (|13|) - (|16D). The symmetry Z (see (^)) does not effect the period T 
of the set, but the symmetry (^) does. Having been applied for odd T = 2?t, + 1, it doubles 
the period: 

P{2/(^)} = {y'-'^'} (30) 

For instance for the period T = 3 

P{yi 2/2 2/3 • • •} = {2/1 -2/3 2/2 -2/1 2/3 -2/2 • • •} 

(We write down only the simplest periodically repeating pattern). 

Note that set {y^^^, odd T for the homogenious model, corresponds to the sequence of 
period 4T for the model with alternating constants: Q{y^^^} = {y^'^^^} 

For the sets with even T = 2n the symmetry (pUj), generally speaking, does not effect the 
period, as can be easily verified. 

To go further, we need an additional piece of information concerning the maximal eigenvalue 
set {y}. We supply it by stating: 

The maximal eigen-value set for the homogenious model is unique, modulo shift, for all 
phases except the ferroelectric phase I (see fig. 3). 

This statement is equal to stating that maximal global transfer matrix eigenvalue is non- 
degenerate. For the 'source' 6- vertex model, this is true: y = for all phases except ferroelectric, 
and y = is just the value invariant under action of Z-symmetry (^). Strong coupling limit 
/i — >• cxD is non-degenerated, too. As to the exclusion - ferroelectric phase I (see fig.3), it has 
degeneracy 2^, K being the number of planes. Indeed the ferroelectric phase is built up from 
a- (or b-) vertexes only. It follows from (|l]) that Boltzmann weight does not depend on h for 
the homogenious model, whatever a- vertexes (type 1 or type 2 -see fig. 2) form each plane. 
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For all other phases, the hypothesis of non-degeneracy means that the maximal set {y}^ax 
is invariant under action of Z— and P— symmetries (^,^), modulo arbitrary shift 
The immediate consequences are: 

1) from the (|25|) : that the averave value of < i/p >max must be zero: 

1 ^ 
^ p=i 

2) from the (|30|) : that the maximal set cannot have an odd period. 

The period of maximal set also cannot be T = 2n, n odd. To see that, consider T = 6. For 
the set to be invariant under action of Z-symmetry (p5D it must have form 

= i^u V w —u —V —w . . .} 

Acting on it by the P-symmetry 

P{y(6)} = ^ -V...}, 

one obtains the set with period 3. Analogously one obtains Pjy^^"^} = {y''"^}, n odd, in 
contradiction with non- degeneracy hypothesis. 

Since the maximal sets {y*-^-*} with T and T/2 odd are forbidden, we are left with the only 
possible choice 

T =4:71 

For that case, 

P{y(4")} = always. 
The first nontrivial case is T = 4. Note that it is just the period which arises in the strong 
interplane interaction limit h oo. According to Z-symmetry (or as well to P-symmetry ) 
invariance, the maximal set reads 

{y^'^ = {v^-v-C--} (31) 

T = 8. P- and Z-symmetry invariance lead to the following set: 

{y(^)} = {u V 0. V -u -V 0. -V . . .} (32) 
With increasing T, admissible structure looks more and more complicated. For T = 12 

= i^u V w r —w V —u —V —w —r w —v . . .} (33) 

and so on. Thus, we have produced the hierarchy of sequences which are candidates for the 
planes polarization set {y}^ax corresponding to the maximal global transfer-matrix eigenvalue. 
The number of parameters is reduced essentially but still there is quite a freedom left. The 
results are in agreement with numerical data. The latter however forces us to formulate the 
final hypothesis which we cannot prove. It is described in the section to follow 



7 The final hierarchy 

Let us assume that the maximal set contains at least one plane k with maximal polarization 
= 1. It follows immediately from (^) and non-degeneracy, that some another plane p has 

^ Note that for the model with alternating constants {...h, —h, h, — /i,...}, the maximal set {y}^cix always degener- 
ate. It is the additional symmetry F that accounts for the degeneracy. The same holds for the model with the constants 
{. . . /i, h, —h, —h, . . .}. 
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Up = —1. Denote the distance between these two planes hj A = \p — k\. We shall parametrize 
the sequences {y} by the value of A. The entries of maximal set {yj^ax obey the following 
rule: 

Vk+n = -Vk^n, n odd 
yk+m = yk-m, rn even 

and the same with the replacement of k by any other number where yj^ = ±1. From (|3^) 
those numbers are k = k i: An, integer n. Making use of (0), for A odd we obtain the y-set 
of period AA {y^^^^}, and {y^"^^^}, for A even. The rules (0) are consistent with Z and P 
symmetries, and can be derived from them. 
Let us take A = 1. Here we have 

{y} = {l 1-1-1...} (35) 

— just the maximal set for the h oo limit (phase 4 in the phase diagram). A = 2 produces 
the set 

M = {i-e-u...} (36) 

- the maximal set for the phase 5 in the phase diagram. It coincides with that from ( ^2]) when 
rj = 1. For A = 3 one obtains from Eqs. (|34D : 

{y} = {1 V w —1 —w V —1 —V —w 1 w —V . . .} (37) 

which coincides with ( p3D when u = 1, r = —1. 

A = 4 produces the set ( |32D with the substitution u = 1, and so on. 

Phases with A = 1,2 do exist on the phase diagram. However we have failed to find 
there numerically the phases with A = 3,4. From that we conclude that another members of 
hierarchy, with A > 4 don't appear, too. We suggest that the higher members of hierarchy do 
appear when one includes into consideration more distant than nearest neighbours, plane-plane 
interaction. For instance, we have checked that the phase with A = 3 (^) does enter the phase 
diagram when one includes additional interaction between the planes k ^ k + 3, any k. 



8 The phase diagram 

We have shown in the last two sections that the set of polarization constants {y} = yi y2 ■ ■ ■ yx 
which corresponds to the maximal eigenvalue (the maximal set) must have the period divisible 
by 4: T = 4n. For the simplest case T = 4 the maximal set is ( pi]) and the transfer-matrix 
eigenvalue is given by Eq. (pOD , with the substitution yi = ys = y^ = rj; y2 = y^ = yK = ^ 



A2/^ = A^(2/iOAe(-2/i?7) 

We can always choose ^ and t] positive. Therefore we only have to maximize the eigenvalue 
([18|) within the two-parameter space 

2n~N 2m -N N 

e = ^^; ^ = ^^; n,m = 0,l...- 

By means of the Newton-Ralphson method we are able to solve BA equations directly and find 
the block with largest eigen-value for system size up to = 32. 
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phas9 diagram of symmetric 6-vertex model, ti=0 



0.5 
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Figure 3: Phase diagram of the 2D 6-vertex model (corresponds to decoupled planes limit h — 0) 



It turns out however that the maximal eigen-value for all values of parameters a, b, c and h 
belongs either to the block n = m = N/2 (in that case h vanishes from all equations (|15|,|16D, 
yielding the known results for 6-vertex model), or to the block with at least one of n, m, say m 
equal to zero: m = 0. Finally the Bethe Ansatz (|15|,|16D rewritten 



A. 



n 



A2/^ = A„=o(2/iOA5(2/i) = A1A2 

_ ^2h(2n~N) _|_ ^N^~2h{2n-N) 

b{2ATj 



-bTj 



n+l TT 



n 



2Ar, - 



-a + br, 



1=1 m - 2AT-i 4 

Note that the eigen-value is invariant under the transformation: 

a b, b ^ a, h ^ —h 



Figures show the phase diagram of the model in the space of Boltzmann weights ratios 
a/c, b/c, for the different values of interplane constant h. Fig. 3 (decoupled planes, h=0) repeats 
well-known results for usual 6-vertex model: the three distinct phases exist separated by the 
lines A = ±1. These phases are (for details see 0): 

I. Antiferroelectric phase. C-Vertexes (type 5 and 6) are dominant; In sufficiently low tem- 
peratures vertex plane configuration is filled with arrows alternating in both directions with 
step 1. 

II. Disordered phase — 1 < A < 1. All types of vertexes are present; there are no types 
of vertexes which are dominant. Here lies the high temperature limit so one should expect 
disorder. 

III. Ferroelectric phase A > 1. Phase III occupies two separated regions, the left upper and 
the right bottom in the phase diagram. In the left upper region, b/c > a/c+1, each plane is 
occupied by either exclusively type 3, or exclusively type 4 vertexes. It is convenient to consider 
the arrow resulting from vector summation of all arrows in the plane: 

for the type 3 configuration, all arrows on the plane are pointing up and to the left, resulting 
is upleft arrow \ , whereas for the type 4 configuration, the resulting is downright arrow \ , 
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phase diagram of SD-extended 6-vertex model, h = 0.15 



0.5 



1.5 2 2.5 3 



Figure 4: Phase diagrams of the SD-extended 6-vertex model on cubic lattice in the a/c, b/c plane, h — 0.15. 
The phase transitions between the phases II/IV, II/V, V/III are of the first order. The phase transitions between 
the phases IV/V are of the second order, Pokrovsky-Talapov type. Transition I/II is the Kosterlitz - Thouless 
type transition. Curves separating phases III/V , V/IV) are given by: A = cosh4/i {h/c = (a/c) exp4/i — 1, 
respectively. 



1, namely, y = 1 for 



see Fig. 2. Polarization vector ([1^) modulus reaches its maximum \y\ 
the type 3 configuration and y = —1 for the type 4 configuration. 

If we associate with the each plane k the corresponding arrow — \ for type 3, \ for 
type 4 configuration, then the state of the multiplane model in phase III, left upper region, is 
characterized by the random sequence 

••■\\\\\\\ \ ••• (38) 

consisting of these arrows randomly placed. Of course switching on some vanishing external field 
favouring type 3 configuration will lift up this degeneracy, giving the homogenious sequence: 



The state of the model in the right bottom region, phase III, a/c > b/c+1 will be analogously 
characterized by the type 1, type 2 vertexes configurations, or upright, downleft y , / 
arrows. 

■■■//////// ■■■ 

In Fig.^ ^ the phase diagrams for the systems with a fixed h in increasing order are given. 
The phases I - III are exactly the same as in uncoupled planes limit in Fig.^. For nonzero 
coupling h the regions occupied by the phase III are given by the formula A = cosh4/i, and 
two new phases arise. 

V. Layered antiferroelectric phase. This phase is exactly the one described as the strong 
interplane coupling limit in section 2. Each plane is in ferroelectric phase, but plane A; + 1 
configuration strictly follows from one in fc-th plane by 7r/2 clockwise rotation. So, the planes- 
layers in increasing order are formed by the unique, up to shifting, sequence of vertexes, or 
resulting arrows: 



(compare with the sequence 

In that phase polarization vector (0) alternates each other plane yk+2 = —yu'- 



■ yk,yk+i, 



(39) 
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Figure 5: Phase diagram of the SD-extended 6-vertex model on cubic lattice in the a/c, h/c plane, h = 0.3 



phase diagram of 3D-extended 6-vertex model, h = D.5 



0.5 1 1.5 2 2.5 



Figure 6: Phase diagram of the 3D-extended 6-vertex model on cubic lattice in the a/c, h/c plane, h = 0.5 
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polarizations distribution in intermediate phase, h=0.15 

' ' «' o»? 'JJ'T • ' 

❖ o 2 <C ^o* o;^mirius1 (x) 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 



Figure 7: Distribution of polarizations in the intermediate phase V. h=0.l5. This phase is characterized by 
the periodically repeated set of polarizations (in planes) ( |3^ ) {1 — ^ — 1 ^ . . .}, ^ = (2n — N)/N, n being the 
number of upward pointing arrows in a row in the corresponding plane. The results of numerical calculations 
of ^ for the system size = 32 are shown. The area is divided into 8 sectors, sector number k corresponds to 
2k ~ 1 < n < 2k. Thus, polarization monotonically increases from ^ = — 1 on the second order phase transitions 
line to ^ = near the point (0,1). On the 'free fermions' circle A = inside the intermediate phase, exact value 
of ^ is given by sm(7r^/2) — —{a/b)sinh{'ih). 



forming the layered antiferroelectric structure (in each separate plane the structure is ferro- 
electric) Polarization vector modulus in each plane reaches its maximal value \yk\ = 1, any 
k. 

V. Phase V is an intermediate phase which can be described as follows: system splits into 
two subsystems, each contains two planes k-th and {k + 2)-th, e.g. (the first and the third) and 
(the second and the fourth). In one of the subsystems the planes are in ferroelectric phase, with 
alternating polarization vector t/k = 1, yk+2 = — 1- Another subsystem is partially disordered 
Uk+i = 2x — 1 = —yk+2,- Again as the previous one, this structure has period 4. Within 
the phase 4, x varies (see Fig. 4a,b); on transition line between phases IV/ V x goes to zero 
continuously, forming the strong interplane coupling structure in phase 4. For the free fermions 
limit, A = 0, or + 6^ = c^, value of x is known exactly (see 0) 

cosTTX = - sinh4/i, for a < b (40) 

Outside the free fermions curve, the branches with equal polarizations are arranged quite 
regular, as is seen from Fig.^,^ However we cannot present the exact formula for the moment. 

The analytical and numerical calculations show that there are no other phases at the phase 
diagram. 

The second order phase transition between the phases V/IV can be found considering the 
equality 

A{n/N ^ 0) = 1 (41) 
. Proceeding analogously as in ^ we get 

. For ae^^ > b,b > a we have 



aTj - b{2ATj - 1) _4h _ T-r a - 26A + be 



4/i 



n ^ r j ^ ^-"=n 



a — bTj ae^" — b 



-e 



-Ah 
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Figure 8: Distribution of polarizations in the intermediate phase V. /i=0.3. 



and I sinh4/i = 1 — Ae or using (p! 



b/c = (a/c) exp (Ah) - 1 (42) 

All phase curves are symmetrical with respect to a = b line. 

Proceeding with the same equality and taking b > ae^^, we obtain: 

A = cosh4/i (43) 

But now this line marks only right margin for the phase transition III/V, and not necessarily 
the exact point. We have found however with all numerical accuracy that line (^31) is exact 
transition point indeed. With increasing h, the phases IV and V expand in space, and the others 
diminish. When h reaches hi = 0.458134 . . . and higher, phase II disappears completely as is 
shown in Fig.^ The fraction of space occupying by the phases II, III, diminishes exponentially 
with h, as is seen from ( ^ , ^3] ). 



Finally, the transition line between the phases II/IV (II/ V) or I/IV (I/V) in Fig.^-|^ is 
obtained numerically. The line ends in points (0,1) and (1,0) which agrees with the limiting 
'Ising chain' case. Phase transition II/IV (II/V) or I/IV (I/V) is of first order. Indeed in phases 
II or I the partition sum does not depend on h at all; (all = N/2 and polarization vector i/k 
for all k is zero. The values of y^-set jump when crossing the critical curve and dependence on 
h shows up. So the partition sum will have a cusp as a function of h. 

The transition IV/V is of the second order. On this line, the order parameters — polar- 
ization vectors i/k ( |14|) change continuously when approaching the critical point. The second 
derivative for free energy over h diverges as inverse square root ~ 1/Vh — h* in the critical 



point (Pokrovsky-Talapov type transition 



We have described the phase diagram of the 3D-extended model with the homogenious set 
of constants. The phase diagrams for the models with arbitrary sets of 'h' and ' — h' constants 
are the same, with redefinition of phases (see Eqs. (p2| - p4D ). Note that for the model with 
alternating constants, {. . .h, —h, h, —h, . . .}, the period is always 2 (in planes). 
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Conclusion 



We have obtained the phase diagram for 3D solvable multilayered 6-vertex model, in full 3- 
parameter space. The model enjoys locality of interactions and positivity of Boltzmann weights. 
The applicability of the method to other solvable vertex models with ice rule ( ^61) is shown in 
0. in view of possible applications, note that the strength of layer- layer interaction h can vary 
from plane to plane, as well as anysotropy parameter within each layer. Another possibility 
is to include more distant than nearest neighbour, interactions along 3-rd axis. The resulting 
solvable models are ones with competing interactions pj. 

Another interesting question is the universality class of the model we have considered. The 



finite size scaling analysis (see e.g. [rT|) of new (due to plane-plane coupling) critical phase, 
named V in the phase diagram, shows that it is described by 2D conformal field theory with 
central charge c = 1. Thus it belongs to the same universality class as the "source" 6-vertex 
model. 
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Appendix 

The partition function for the system with open boundaries is given by 

Z = T'' 

where T is the global monodromy matrix. In our case (see Eqs. (^ |TT])) T factorizes into 
product 

'r=f[% (44) 

k=l 

where the 'local monodromy matrix' for the fc-th plane 

n=l 

Ln being the matrix of Boltzmann weights for the 6-vertex model, upper index (k) corresponds 
to the k-th plane, matrices a and r are defined by (|12]). Omitting the index k in the right-hand 
side of the last formula for convenience and denoting 

^ ^. r^k+i) (fc-i) / // ^ hi. ^ -^a 
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we have: 

N 



n=l 



-hcTT- 



1 Li e^"< e-^""2 L2 e3""2 . . . g"'^--^ Ln e^"^'^ (45) 

The exponential factors in this expression commute with each other because they are diagonal 
matrices (see (0)). The commutation with L„ is given by 

[exp {h{a + Tn)Tp)), Ln] = [exp (/i(cr + r„)rp)), L„] = 0; 

which is equivalent |^ to the charge conservation property of the "source" 6- vertex model 
i): 

L^va^R' = 0, unless a + (3' = a + (3 (46) 



Using these commutation rules, we move all exponents in (^51) outside to the left and to the 



right. For instance, to move the term e ^"'^^ to the left, one inputs unity 



The first factor in this expression commutes with Li and goes to the left while the second 
one commutes with all Li, L2,. . .L^ (because they act in the different subspaces) and goes 
to the right. Repeating the similar procedure for all exponents, and taking into account 

N m~l N-1 N 

E E = E E ' we obtain 

m=2 n=l n=l m=n+l 

%{Kg) = e-^'^ Sn^i^"' T,(0,0) e^-^ Sli--" A~^^ (47) 

( N m~l \ 



\m=2 n=l 

N 



Then, due to the charge conservation (^61) , ^ r^'^-' is a constant. For the 6- vertex model 



n=l 



this property is known also as the ice rule (see e.g. |]): 



1 ^ 

n=\ 



the so-called polarization vector. may take values — 1 < < 1- 
Restoring the index /c, can be written as 



N m-1 



^m=2 n=l 



Multiplying 7^7^ ■ ■ - Tk and using A^+i = C*^ Ck+2,k+i, we get for the global mon- 
odromy matrix (^: 

r = Ci,oC2,i e-'^^'-'''^' ri(0, 0) e'^o'^^''^" C^iC;,2 . . . 
Ck-i,k-2 Ck+i,k e-'^-'^'"'^-- r^(0,0) e'^--'^'"'^-- C^^,^^ C],]^_, 
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Here Yj, = N y^. 

Moving C-factors (C~^-factors) to the left (to the right), one obtains 



K 



k=l 

B depends only on r^'^^ and therefore it is a gauge transformation. In the other hand, one 
can consider the transfer-matrix for the 6- vertex model in a horizontal field of strength H 



N 



n=l 



Doing the same procedure, we arrive to formula (|47|) with substitution h = —g = H/2; r^, r^' 
1: 

in which A again is a gauge transformation. Comparing the last relation with (0), we see that 
partition function of our model and the set 6-vertex planes, fc-th plane in a field 

Hk = hkVk+i — hk-iVk-i 

(both systems with open boundaries) are gauge-equivalent. Hence, their transfer-matrices are 
equal 

K 

^ our model 

n TfSvihkVk+i - hk^iVk-i) 

k=l 
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